Stability and eccentricity of periodic orbits for two planets in a 

1:1 resonance. 
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ABSTRACT 

The nonlinear stability domain of Lagrange's celebrated 1772 solution of a three-body problem 
is obtained numerically as a function of the masses of the bodies and the common eccentricity 
of their Keplerian orbits. This domain shows that this solution may be realized in extra-solar 
planetary systems similar to those that have been discovered recently with two Jupiter-size planets 
orbiting a solar-size star. For an exact 1:1 resonance, the Doppler shift variation in the emitted 
light would be the same as for stars which have only a single planetary companion. But it is more 
likely that in actual extra-solar planetary systems there are deviations from such a resonance, 
raising the interesting prospect that Lagrange's solution can be identified by an analysis of the 
observations. The existence of another stable 1:1 resonance solution which would have a more 
unambiguous Doppler shift signature is also discussed. 

Subject headings: extra-solar planetary systems, resonances, Lagrange solution 



1. Introduction 

In many extra-solar planetary systems discov- 
ered recently, the observed Doppler shift of the 
emitted light is well described by assuming that 
the central star is moving in a Keplerian elliptic 
orbit due to its gravitational interaction with a sin- 
gle Jupiter-size planet. If residuals are present in 
a least squares fit to the data after possible chro- 
mospheric fluctuations in the star have been taken 
into account, these residuals signal the presence of 
an additional planet, or possibly several planets 
(Fischer et al 2000) (Marcy et al 2001a). There 
is, however, an important exception when such a 
star also travels on a Keplerian orbit even though 
there are two planets orbiting it. We have in mind 
Lagrange's celebrated solution of the three-body 
problem (Lagrange 1772), for which he won the 
prize of the Royal Academy of Science of Paris. 
In this solution, each body moves on a Keple- 
rian orbit with a common plane, period, eccen- 
tricity and focus that is located at their center of 
mass, in such a manner that at all times the rel- 
ative positions of these bodies form the vertices 
of an equilateral triangle of variable size, (see Fig. 



1). This solution will be discussed in section 2. 
It is therefore interesting to consider the possible 
occurrence of such a 1:1 resonance in extra-solar 
planetary systemm with two Jupiter-size planets, 
particularly in view of the recent discovery of a 
remarkable 2:1 resonance of two such large plan- 
ets in GJ876 (Marcy et. al. 2001b). To be rele- 
vant to current astronomical discoveries, however, 
it is necessary that Lagrange's solution be stable 
in the range of observed masses and eccentrici- 
ties. In the past, linear stability analyses have 
been carried out that were primarily focused on 
the restricted three-body problem (where one of 
the three masses vanishes), which is applicable to 
the study of the motion of asteroids, but these re- 
sults can also be extended to the general Lagrange 
solution as will be discussed in section 3. In this 
paper we present the results of a numerical analy- 
sis of the nonlinear stability domain of Lagrange's 
solution as a function of the masses of the three 
bodies and the eccentricity of the common ellipti- 
cal orbits (see section 4). In practice there will be 
deviations from Lagrange's solution, and the inter- 
esting question arises as to whether there are dis- 
tinct characteristics which would distinguish these 
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deviations from other types of perturbations due 
to additional planets. While variations in the ec- 
centricity and major axis of the approximate ellip- 
tical orbit are common to all perturbations caused 
by the presence of a second planet, one of the most 
distinguishing feature of a slightly off-resonance 
but stable Lagrange solution is that the the rota- 
tion rate of the axis of the elliptic orbit of the star 
is much smaller than for other types of perturba- 
tions, (see Fig. 2). 

For example, for a 1% deviation in the posi- 
tion of the lighter planet, assumed to be about a 
Jupiter mass, we find that it takes about 800 pe- 
riods for a complete revolution of the major axis, 
(see Fig. 2). In contrast, a fit to the recently dis- 
covered 2:1 resonance in GJ876 indicates that the 
major axis of the planets should complete a rev- 
olution in about 53 periods of the heavier planet 
(Nauenberg 2002). 

Lagrange's solution is not the only 1:1 reso- 
nance for the three-body problem. We also con- 
sider another solution that we found to be stable 
in the domain of masses and eccentricities relevant 
to extra-solar planetary systems (Laughlin 2002 ) . 
In this case the two lighter bodies (planets), and 
the heaviest body (star) have different orbits, and 
therefore the Doppler shift data should be readily 
distinguishable from the case of a single planet. 
The characteristic feature of this solution is that at 
each half period the three bodies are aligned, and if 
the heavier planet is in a nearly circular orbit, the 
lighter planet moves in a highly eccentric elliptical 
orbit, (see Fig. 3). A periodic alignment of the two 
planets and the star is also a characteristic of other 
resonance solutions, as in the case of the 2:1 res- 
onance in GJ876 (Laughlin and Chambers 2001) 
(Rivera and Lissauer 2001) (Lee and Peale 2002) 
(Nauenberg 2002). The evolution of this configu- 
ration when the planets are slightly off resonance 
is shown in Fig. 4, which illustrates the clockwise 
rotation of the major axis of the eccentric orbit. 

2. Lagrange's solution of a three-body 
problem 

In Lagrange's solution (Lagrange 1772) of a 
three-body problem, each body travels on a sep- 
arate elliptic orbit with a common period, eccen- 
tricity and focus which is located at their center 
of mass, in such a way that these bodies are al- 
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Fig. 1. — Lagrange's periodic 1:1 resonance solu- 
tion for the three-body problem, showing the loca- 
tions of the bodies at apogee on the vertices of an 
equilateral triangle (full lines). The dashed lines 
show the direction of the major axis for each el- 
lipse, and their intersection at the center of mass, 
which is their common focus. 
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Fig. 2. — Rotation of the axes of the ellipses in 
Lagrange's solution after 800 periods, for a one 
percent deviation in the initial position of planet 
2 from exact resonance. 
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Fig. 3. — Another 1:1 resonance solution for 
the three-body problem, showing the inner heav- 
ier planet (green) on a nearly circular orbit, the 
lighter planet (red) on an eccentric elliptic or- 
bit with eccentricity e = .8, and the central star 
(blue), as they appear aligned at maximum elon- 
gation of the ellipse. 




Fig. 4. — Evolution of the orbits shown in Fig. 
2 for a small deviation from 1:1 resonance, show- 
ing the rotation of the major axes of the planets 
during 10 periods. 



ways at the vertices of an equilateral triangle of 
variable size. An example is illustrated in Fig.l, 
where the masses are in the ratio 1.0 : 0.7 : 0.2, 
and the common eccentricity is e = .5, showing the 
equilateral triangle for the relative positions of the 
three bodies at apogee. Apart from the overall ori- 
entation of the system, these parameters uniquely 
describe the solution. In this solution the rela- 
tion between variables is somewhat different from 
the corresponding ones in the two-body problem. 
For example, the common frequency to or period 
P = 2tt/lu of the motion is given by the relation 



w= y/Gm(l + e) 3 /R 3 , 



(1) 



where G is Newton's constant, m — m\ -\-m 2 + m 3 
is the total mass, mi is the mass of the heaviest 
body (star), m 2 , m 3 are the masses of the lighter 
bodies (planets), e is the common eccentricity, and 
R is the maximum size of the equilateral triangle 
on which the three bodies are located. Then the 
major axis ai of each of the elliptic orbits is 



+ + rrijiJik 



m(l 



-R, 



(2) 



where the subscripts i,j and k are permutations 
of the the integers 1,2,3, while the mean of the 
velocity of the star at the maximum and minimum 
distance from the foci of the ellipse is 



K = (2irG/P) 



1/3 



V m 2 + 



m 2 m 3 



2 / 3 Vl 



(3) 



If the planetary masses TO2 and 7713 are small com- 
pared to the mass of the star, ml, then a 2 ~ 03 w 
a, where a = Rj (1 + e) and uo w \J Cm/ a 3 as in 
the corresponding two-body problem. In princi- 
ple, any fit to the data with a single planet of mass 
vrip can also be attributed to two planets which, 
according to Eq. 3 for K, have masses satisfying 



the relation 



= Vr 



m 2 m 3 . 



which is 



somewhat less than the sum of the masses of the 
two planets. In practice, however, it is unlikely 
that extra-solar planetary systems would occur in 
an exact 1:1 resonance, and therefore the presence 
of a second planet manifests itself in the occur- 
rence of residuals in a single Keplerian orbit fit to 
the data. 



3 



3. Linear stability analysis 

The first linear stability analysis of Lagrange's 
solution was carried out by Routh for the special 
case of circular orbits (Routh 1875) 1 . Assuming 
that the attractive forces between the bodies de- 
pends on the relative distance r as l/r K , Routh 
demonstrated that Lagrange's solution was stable 
provided that the masses satisfied the inequality 



7< 3 ( TT^' 



(4) 



where 



(mim 2 + mim 3 + m 2 m 3 ) ,_. 
(mi + m 2 + m 3 ) 2 

For gravitational interactions where k = 2, the 
constant on the right hand side 2 of Eq. 4 is 
1/27, which implies that the masses of the two 
lighter bodies (planets), m 2 and 7713, must be much 
smaller than the mass of the heaviest body (star) 
mi. For example, setting m 2 = 0, which corre- 
spond to the restricted three-body problem and 
applies to the motion of asteroids such as the Tro- 
jans, this inequality implies that 777,3/(7711 +1713) < 
.03852.. which has become known as Routh's crit- 
ical point. Neglecting quadratic terms in the mass 
ratios m 2 /mi and m 3 /mi, Routh's inequality be- 
comes approximately (m 2 + m 3 )/mi < 1/27, in- 
dicating that the stability depends to a very good 
approximation only on the sum of the masses of 
the lighter bodies. 

For the stable configurations, Routh obtained 
the frequencies uj\ , u 2 and uj 3 of the normal modes 
in the plane of the orbits, which for k = 2 are 
given by 



^(1 + v/I - !^, 



LJ 3 = U), 



(6) 



(7) 



(8) 



1 In this paper, Routh incorrectly attributed Lagrange's so- 
lution to Laplace. 

2 Routh found that he had been anticipated in this impor- 
tant result by someone called M. Gaschcau in a thesis on 
mechanics, apparently left unpublished 



where u> is the fundamental Kepler frequency Eq. 
1, and he determined the corresponding ampli- 
tudes for these modes. While Routh did not dis- 
cussed the stability with respect to deviations per- 
pendicular to the plane of the orbit, it is straight- 
forward to show that for this case the Lagrange 
orbits are stable. Remarkably , he also briefly 
considered second order deviations, which he re- 
marked could "ultimately disturb the stability", 
but this part of his analysis was incomplete, call- 
ing attention only to two of the possible commen- 
surability or resonance relations, 0J3 = 2w 2 and 
u)\ = 2u>2, for which 7 = 1/36 and 7 = 16/675, 
respectively. 

For elliptic orbits, a linear stability analysis 
of Lagrange's solution was not carried out un- 
til some 90 years later, when (Danby 1964) nu- 
merically integrated the Floquet equations for the 
first order deviations from the Lagrange solution 
of the restricted three - body problem (see also 
(Bennett 1965). Subsequently, a majority of sta- 
bility studies have been confined to this special 
case, but as we shall see, the results can also be 
applied to the general solution (Marchal 1990)). 
For a modern discussion of the stability of the La- 
grange solution see (Siegel and Moser 1971) 

4. Nonlinear stability domain of Laplace's 
solution 

The nonlinear stability domain of Lagrange's 
solution for the three-body problem that is pre- 
sented here was obtained by integrating the equa- 
tions of motion numerically, and determining 
whether, for a small initial deviation from the 
solution, the orbits were either confined or uncon- 
fincd after a large number 77 of periods. We found 
that when 77 was increased from 400 to 800, there 
were no significant changes in the results except 
near the critical points discussed below, where we 
increased 77 until no further changes occurred. For 
the deviations we consider small displacements of 
the position of one of the lighter bodies planets 
in the plane of the orbit and also perpendicular 
to this plane. As we shall see, except in the case 
when o>3 = 2cj 2 , only initial deviations in the plane 
of the orbit gave evidence for some of the expected 
nonlinear instabilities near the commensurability 
relations for the three fundamental frequencies, 
Eqs. 6-8. 
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Starting with a small deviation in the velocity 
of one of the lighter bodies of order 1% in the direc- 
tion perpendicular to the plane of the orbits, our 
numerical result for the nonlinear stability domain 
is shown in Fig. 5, where the unstable regions are 
indicated by small squares. In agreement with the 
linear stability analysis ( see section 3) , the nonlin- 
ear stability domain depends only on the common 
eccentricity e of the orbits (horizontal axis), and 
the Routh parameter 7, Eq. 5 (vertical axis). This 
parameter is determined by the ratios m 2 /mi and 
7713/7711 of the masses of the lighter bodies (plan- 
ets) to the mass mi of the heaviest body (star) 
(vertical axis). We have verified this result by 
computing this domain for different fixed values 
for the ratio of the masses of the planets, without 
observing any changes when plotting the results 
with Routh's parameter 7. For the particular case 
shown in Fig. 5, the computation was carried out 
for m 2 = 7713, while the curves shown are the lin- 
ear stability computations of Danby and Bennett 
(Danby 1964) (Bennett 1965), which were origi- 
nally obtained for the restricted three-body prob- 
lem where either 7772 or 7773 is set equal to zero. 

This nonlinear stability domain looks surpris- 
ingly similar to the stability domain of the Math- 
ieu equation (Grimshaw 1990) with an additional 
nonlinear restoring or damping term. Indeed, for 
the restricted three-body problem, a linear stabil- 
ity analysis of the equations of motion, to first 
order in the common eccentricity e, in a frame of 
reference rotating with the frequency of the or- 
bits, yields a bifurcation at u 3 = 2u 2 which can be 
viewed as a parametric resonance between the fun- 
damental orbital frequency for elliptic motion and 
the oscillation frequency of the first order devia- 
tions of the massless body. In the general case, we 
have seen that this bifurcation occurs for e = at 
7 = 1/36 = .02777.... For the restricted problem, 
77i2 = 0, this corresponds to A = 2\/2/3 w .9428, 
where A = (mi — m^)/{mi + 7773) and 7773/7771 w 
.02944. We evaluated the instability domain to 
first order in the eccentricity e, and found that it 
lies inside the wedge 

A - Se < A < A + 5e, (9) 

where <5 = 3/20\/2 w .10606. This domain is 
shown in Fig. 6, together with the correspond- 
ing nonlinear domain where the orbit actually be- 
comes unstable. These linear boundaries mark 



the onset of a bifurcation in which the orbits first 
become aperiodic filling out a confined region of 
space. 

An enlargement of the stability domain near the 
critical point at e = and 7 = 1/27 w .037037.. 
is shown in Fig. 7, which displays our results 
for the case that m 2 = 7773. This critical point 
is often discussed in connection with the stabil- 
ity of the Lagrangian points for the restricted 
three body problem in a frame of reference rotat- 
ing with the frequency of the circular orbits, and 
is applied to the study of the Trojan and other 
asteroids (Murray & Dermott 1999). As in the 
previous case, however, this critical point marks 
only a bifurcation to an aperiodic but confined 
motion, while numerically we find that the non- 
linear instability begins at a somewhat higher 
value of 7 w .0391, corresponding for m 2 = to 
777.3/7771 w .0425. A curve quadratic in the eccen- 
tricity fits well the boundary of the upper part 
of the nonlinear instability domain, as shown in 
Fig. 7, while a linear curve fits the lower part of 
the domain with a slope which is somewhat higher 
than the one which we calculated analytically in 
the linear approximation, Eq. 9 but in good agree- 
ment with the numerical linear stability results of 
(Danby 1964) and (Bennett 1965). 

Up to now, we have considered only initial devi- 
ations in the direction perpendicular to the plane 
of Lagrange's orbits. In Fig. 8 we show the nonlin- 
ear stability domain obtained by taking an initial 
a displacement dx = dy = .001 in this plane for 
one of the two lighter bodies. This result was ob- 
tained by fixing the ratio 7772/7773 = 1, but the 
same results are obtained for other values of the 
ratios of these masses when the results are plot- 
ted as a function of the Routh parameter 7, Eq. 
5. We see evidence in this figure for the insta- 
bilities at the nonlinear resonances u)\ = 2cj 2 and 
lu 3 = 3ut2, corresponding to 7 = 16/675 w .02370.. 
and 7 = 32/2187 w .01463.. Also, the domain 
of stability shrinks in the upper wedge with in- 
dications of additional resonances in this region. 
This is confirmed by an enlargement of this re- 
gion, which shows the tail of the instability due 
to the resonance at 3u>i = 4w 2 corresponding to 
7 = 576/16875 w .03413... as well as other reso- 
nances (see Fig. 9). 
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Fig. 5. — Stability domain for Lagrange's solution 
as a function of the eccentricity e and the Routh 
parameter 7, Eq. 5, obtained for initial deviations 
perpendicular to the plane of the orbit. 



Fig. 7. — Enlargement of the stability domain 
shown in Fig. 5 near the critical point at u)\ = 002- 
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Fig. 8. — Stability domain for Lagrange's solution 
Fig. 6.— Enlargement of the stability domain as a function of the eccentricity e and the Routh 

shown in Fig. 5 near the critical point for the parameter 7, Eq. 5, obtained for initial deviations 

parametric resonance at uj 3 = 2lu 2 - 111 tne P lane of tne orb 1 *- 
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5. Another 1:1 resonance solution 



Lagrange's solution is not the only stable 1:1 
resonance system for the three-body problem. An- 
other solution is illustrated in Fig. 3, in which the 
two lighter bodies or planets are traveling on two 
different orbits which are approximately elliptic in 
such a way that the central star and the planets 
are aligned when located at the maximum or min- 
imum distance from the center of mass. We found 
these periodic orbits (Laughlin 2002 ) by a new 
method based on an expansion of the coordinates 
in a Fourier series with a common period, and we 
determined the Fourier coefficients of this expan- 
sion by finding minima of the action integral with 
respect to these coefficients by an iterative pro- 
cess (Nauenberg 2001). A characteristic feature 
of these orbits is that when the heavier of the two 
planets is in a nearly circular orbit, the lighter 
planet is in a elliptical orbit with a very large ec- 
centricity which rises slowly with increasing mass 
of the heavier planet, as shown in Fig. 10 . 

In contrast to the Lagrange's solution discussed 
previously, a small deviation from exact 1:1 reso- 
nance in this case leads to a relatively rapid ro- 
tation of the major axis of the elliptical orbit, as 
illustrated in Fig 3. This leads to characteristic 
modulations in the Doppler shift oscillations of 
the light emitted by the star. As in the case of 
Lagrange's solution discussed in the previous sec- 
tion, we found that these orbits are stable over the 
range of masses relevant to extra-solar planetary 
systems, see Fig. 11. 

6. Conclusions 

The nonlinear stability domain for Lagrange's 
solution of the three body problem shown in Fig. 8 
indicates that there is a wide range of Jupiter-size 
planetary masses (including brown dwarfs) and ec- 
centricities for which such solutions can exist in 
extra-solar planetary systems. For example, for 
an eccentricity of e w .6 the ratio of the total 
mass of the two planets to the mass of the star 
for which the solutions are stable is .004, except 
for a small region where nonlinear resonances oc- 
cur. This mass correspond to 4.2 Jupiter-masses 
for a one solar-mass star, while for smaller eccen- 
tricities, e < .2, there is a wedge of stable solu- 
tions for higher mass ratios up to approximately 
.04. In principle, any Doppler shift data that can 
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Fig. 9. — Enlargement of the stability domain for 
Lagrange's solution in the upper wedge shown in 
Fig. 8 as a function of the eccentricity e and the 
Routh parameter 7, Eq. 5, obtained for initial 
deviations in the plane of the orbit. 
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Fig. 10. — Dependence of the eccentricity of the 
lighter planet on the mass of the heavier planet for 
the 1:1 resonance solution shown in Fig. 3. 
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be fitted under the assumption that only a single 
planet is orbiting the central star, can equally well 
be attributed to two planets orbiting the star ac- 
cording to Lagrange's solution. In practice, how- 
ever, it is very unlikely that two planets are in an 
exact 1:1 resonance, and therefore one expects to 
find residuals in the data which signal the present 
of a second planet. One of the main effects due 
to the perturbations caused by a second planet is 
the secular rotation of the major axes of the ap- 
proximate Keplcrian ellipses that characterize the 
orbits of the planets and the star. This effect, 
however, would not occur in the case of an ex- 
act 1:1 Lagrangian resonance, and it is strongly 
suppressed in the case that the resonance is ap- 
proximate. It may be thought that in view of the 
greater number of degrees of freedom present with 
Lagrange's solution, a better least squares fit to 
the data should be readily available. This is not 
the case, because Lagrange's equations also allow 
for unphysical solutions where the mass of one or 
even both of the planets can have negative values, 
provided the total mass is positive. Indeed, in a 
preliminary attempt to obtain a least squares fit 
to data which show residuals, the optimal mass of 
the smaller of two Lagrange planets turned out to 
be negative, which ruled out this solution. 

The stability domain for the type of 1:1 reso- 
nance solution discussed in section 5 is shown in 
Fig. 11, demonstrating that this solution also en- 
compasses the possibility of two Jupiter-size plan- 
ets orbiting a solar-size star. In this analysis we 
restricted the heavier planet to be in a nearly cir- 
cular orbit, and found that the lighter planet is 
in a highly eccentric orbit with e w .8 , (see Fig. 
4). If this restriction is relaxed, we also find simi- 
lar stable solutions, and for equal masses the two 
planets exchange eccentricity when the major axis 
rotates through 180 degrees. 

In summary, it is likely that extra-solar plane- 
tary systems that have several Jupiter-size plan- 
ets which are close enough to give rise to signif- 
icant gravitational perturbations will be in res- 
onance, because numerical investigations have 
shown that such systems can be stable. In such 
cases, the planets and the central star are period- 
ically aligned. An interesting exception is the 1:1 
resonance solution of Lagrange, where the planets 
and the star are located at all times on the vertices 
of an equilateral triangle of varying size. It would 



be very exciting if this solution, discovered by La- 
grange 230 years ago, and realized thus far only 
in the motion of the Trojan and other asteroids 
in our solar system (Murray & Dermott 1999), 
would be also present in the orbits of planets in 
extra-solar systems. Likewise a search should be 
undertaken to find also two planets in extra-solar 
systems which are in a 1:1 resonance of the type 
discussed in section 5, which does not occur in our 
system. 

I would like to thank Richard Montgomery and 
Carles Simo for useful references to the vast lit- 
erature on the Lagrange solution of a three-body 
problem. 

michael@mike.ucsc.edu. 
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